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Abstract 


A theoretical model based on laminar boundary layer 
flow equations is developed to predict the erosion rate 
of a graphite (AGCarb-101) surface exposed to a hot super- 
sonic stream of hydrogen gas. The supersonic flow in the 
nozzle outside the boundary layer formed over the surface 
of the specimen is determined by assuming one-dimensional 
isentropic conditions, An overall surface reaction rate 
expression based on the experimental studies by Clajrke 
and Fox is used to describe the interaction of hydrogen 
with graphite. A satisfactory agreement is found between 
the results of the computation, and the available experi- 
mental data. Some shortcomings of the model, and further 
possible improvements are discussed. 
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Introduction 


In order to achieve performance beyond that attain- 
able by the best chemical rochet., the nuclear reactor has 
been applied to a rocket propulsion system. Over the past 
decade, considerable progress has been acheived in this 
direction, and the so-called NERVA (Nuclear Engine for 
Rocket Vehicle Application) propulsion system has received 
considerable development effort. Currently the feasibility 
of a similar, but smaller nuclear rocket engine is being 
explored. Thp NERVA flight engine has a nozzle extension 
after the regenerat ively cooled nozzle,^" and the nozzle ex- 
tension is made of a graphite composite (AGCarb-101) mater- 
ial having the necessary structural strength as well as low 
secondary radiation qualities to reduce shielding require- 
ments. A subscale experimental study was conducted in order 

to obtain the erosion rates of graphite composites when ex- 

2 

posed to supersonic hydrogen atmosphere. We have developed 
a theoretical model based on laminar boundary layer flow 
equations to predict these erosion rates. 

Some details of the experimental setup can be found 
in Reference 2. The continuous supply of high pressure, 
hot hydrogen gas was obtained by using a one megawatt 
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plasma arc generator. The NERVA engine conditions were 
simulated by letting the gas expand through a specially 
designed water-cooled conical nozzle, and placing the 
specimen at a fixed position downstream of the throat. 

The duration for each test varied from half an hour to 
one hour. The heat and mass transfer as well as the chemi- 
cal reaction take place within the boundary layer region 
over the surface of the specimen, and result in the erosion 

of the specimen. Due to the small magnitude of the ob- 

2 

served erosion rates, it is reasonable to treat the bound- 
ary layer flow as steady. Furthermore, for preliminary 
analysis, the flow will be considered laminar, and the 
presence of any dissociated component of the gas in the 
chamber will not be included. The kinetic data dealing 
with the reaction of a graphite composite (AGCarb-101) 
with hydrogen is not available. Therefore, the measure- 
ments of the reaction rate of a graphite filament with 

3 

hydrogen by Clarke and Fox are incorporated in this 
model. It was concluded by these authors that the re- 
action is a surface reaction below 3000°K. The flow out- 
side the boundary layer is determined by assuming one - 
dimensional isentropic conditions. The details are des- 
cribed in the following sections. 

Governing Equations 

Assuming Fick's law of binary diffusion, the steady- 
state axi-symmetric laminar boundary layer flow for a 
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gaseous mixture is described by the following set of 

. ■ 4 

equations 

Mass conservation: 


5 (pur) a (pvr ) 
Bx By 

Momentum conservation: 


( 1 ) 


pu 



+ pv 




Energy conservation: 



Species conservation: 


pu 



+ pv 



JL 

ay 



dy ) 


( 2 ) 


( 3 ) 


( 4 ) 


In addition, we use the ideal gas equation of state for 
the mixture, that is, 

P = p r°t (e h-) 

The pressure distribution in the boundary layer is deter- 
mined from the non-viscous flow outside the boundary layer, 
and under the assumption of one-dimensional flow in the 
nozzle, is given by 



= p u 
p e e 


/ du e \ 
\dx/ 


( 6 ) 
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We are thus dealing here with essentially five unknowns; 
namely, p, u, v, T, and Y^, satisfying five equations 
(Eqs. (1) to (5)). A schematic representation of the dis- 
tance coordinates (x,y) and r as well as the velocity com- 
ponents (u , v) is shown in Figure 1. 

The velocity, temperature, and species mass fractions 
are subject to the following boundary conditions. 

At y = 0 (eroding surface) : 




u = 

0 


(7) 



V = 

v = 
w 
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(8) 

k m 

w \dy/ 

J w 

- E 

i 

h. 

1W 

f 

uu . 

1W 

t 

■ [ k s (W)] - 

2 w 

(9) 

0) . = - 
1W 

f 


D (■ 
w \ 

™±) 

ay'w 

+ UL) Y. 

cw 1W 

t t 

(10 


At y = * (boundary layer edge) : 

u = u (11) 

e 

T = T (12) 

e 

Y r =1, and Y, = 0, i ^ H 2 (13) 

The subscript w ” refers to the conditions interior of the 
specimen (s) . The initial conditions as well as the de- 
tailed expressions for the transport parameters (q, c ., 

li 1 

k, D) are discussed later. 

3 

From the experimental work of Clarke and Fox, it is 
reasonable to conclude that the dominant product in the sur- 
face reaction of hydrogen with graphite in the temperature 
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and pressure range of interest is the methane gas. Then, 

it is necessary to solve only one species conservation 

equation, say for Y„ , since the following identity holds 

2 

good in the mixture, 




1 


(14) 


Since it is convenient to deal with the 'locally simi- 
lar' form of the boundary layer equations, we introduce 
new independent variables by combining the Levy and Mangier 
transformations and the Howarth-Dorodnitzyn transformation 
Let 


s 


.o (p e lI e ll e r2 > dx ' 


[gm 2 sec -2 ] 


(15) 


and 



[dimensionless] (16) 


We also define non-dimensional x-component of velocity (U) , 
temperature (0), and a modified transverse velocity (V) as 
given below 


U = (u/u e ) 
0 = (T/T 0 ) 


V = 


/2 s 


p u u r 
H e h e e 


<-> + (^PXtS) » 


In terms of the new variables defined by Eqs. 
to (19) , the governing Eqs. (1) to (4) become 

+ (e) 


(15) 


(17) 

(18) 

(19) 


( 20 ) 
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The boundary conditions stated in Eqs. (7) to (13) trans- 
form into the following forms: 

At ri = 0 (eroding surface) : 


U = 0 



At ri = 00 (edge of the boundary layer) : 

U = 1 
6 = 1 

Y h =1, and Y i = 0, i ^ H 2 


(21) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


(30) 



Note that the energy boundary condition at the eroding 
surface requires a solution of the heat transfer problem 
within the specimen. Except for the unexposed end por- 
tions, where the heat loss became important, the back sur- 
face of the specimen was heated by an ATJ graphite heating 

element in the experiments conducted for the determination 

2 

of the erosion rates at Aerojet Nuclear Systems Company. 

A realistic treatment of heat transfer within the specimen 
for this situation is too complex, and will also involve 
some unknown parameters. It is, therefore, proposed to 
assume a known constant value of the surface temperature 
(9 ) which replaces the boundary condition stated in Eq. (26). 
Thus , 

9 (s , 0) = e w (3 

The agreement between the theoretical and the experimental 
results will be improved by making different choices of 9 . 
External Flow 

The flow outside the boundary layer over the specimen 
is dependent on the chamber conditions as well as the nozzle 
characteristics. The slight amount of dissociation present 
at the chamber conditions will be neglected. It is reason- 
able to treat the flow in the slender nozzle as one dimen- 
sional. Since there occurs a significant change in temper- 
ature along the nozzle, the variation in Y due to tempera- 
ture cannot be completely ignored. To simplify the calcu- 
lations, we will use two constant average values of Y such 
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that one of them (y-j-) is applicable for the region starting 
from the chamber to the leading edge of the specimen while 
the other (y ) is for the flow over the specimen. 

Therefore, for steady, one-dimensional and isentropic 
conditions, the velocity field in the flow external to the 
boundary layer over the specimen surface can be determined 


from the following differential equation, 


2 u sin a • 
e 


f)'l 


7+1 _ /YPg 

V 


'P Y 
K o 


(p u r 2 1 1 

l K o o o 

L (r Q + x sin a) 2T_1 "* 

' (PnUo r o 3 ) V ~^ 1 

■ (r + x sin a) 2 ^ 2 
v o 


In terms of the boundary layer variable s, defined in Eq, 
(15), we can show that 


P (s) 


p u u L 
K e H e e 


due \ 
dx j 


r(s) = -I 


(Y-l) u. 




PoN ^ 1 


P (s) 


The density (p g ) in the external stream will be obtained 
by making use of the mass conservation, that is. 


p = (p u r 2 )/(u r s ) 
^e H o o o e 


The pressure (p ) , and the temperature (T ) can be cal- 
culated from the isentropic conditions given below. 


P\ /P\Y ,T s \Y-1 

£ -(r) -(£) 
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Note that 

r = r Q + x sin a = r Q + x L sin a (37) 

The only other variable of interest in the external stream 

is the coefficient of viscosity for pure hydrogen gas (d e ) > 

[see Eq. (15)], which will be evaluated as a function of 

temperature (T ) as described in the next section, 
e 

Transport Properties 

The temperature dependences of the specific heats at 

constant pressure for hydrogen (c „ ) , and methane (c ) 

P' ti 2 P* cH 4 

gases were described by the expressions given in McBride et al.^ 

The coefficients of viscosity, and thermal conductivity for 

hydrogen and methane gases were calculated from the expressions 

based on kinetic theory, and the properties of the mixture 

were then determined by making use of the empirical equations. 

Similarly, the binary diffusion coefficient (D) was calculated 

by using the kinetic theory expression quoted in the previously 

cited reference. The effect of the collision integral param- 

(1 1 ) * 

eters appearing in the various transport coefficients, 0 ^ 

(2 2 ) * 

and fK ' , was also included by approximating their temper- 

ature dependences by different analytical expressions valid 
over suitable ranges of temperature. The empirical expressions 
obtained in this manner were always within 2-3% of the values 
obtained from the "true" curves. It must, however, be noted 
that, due to the very small rate of erosion observed in the 
experiments, the mixture will consist of mostly hydrogen gas. 


$ 
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and, therefore, the transport coefficients will be insensi- 
tive to the mixture composition but will nevertheless be 
functions of temperature. 

Chemical Kinetics 

3 

Clarke and Fox concluded on the basis of their experi- 
mental studies that below 3000° K the reaction of graphite 
and hydrogen between 0.01 and 1.0 atm was a surface reaction 
whose rate was proportional to the hydrogen pressure and 
the square root of the dissociation constant of hydrogen. 

An empirical equation for the reaction rate, obtained by them 
and found to be in reasonable agreement with observed values, 
can be expressed as follows: 



-4 -2 

where k c = l - 20 x 10 moles of carbon cm sec 
and kp, the dissociation constant, is given by 


(atms) 


-3/2 


k D A D exp ( E d /R T w ) 


(39) 


with A d = 1*947 x 10 6 atm, and E D = 1*086 x 10 5 cal mole "K 
Note that Clarke and Fox assume in their analysis of the 
experimental data that chemical equilibrium exists at the 
surface of the specimen. In view of the assumption of 
a single step overall reaction, 

C (s) + 2H 2 - CH 4 


(40) 



11 


we have the following relationship 


uu 


h 2 


= (2W„ 

, w H, 


/W ) 


0 ) 


c ,W 


(41) 


The mathematical formulation of the model is now com- 
plete, and a numerical solution is obtained as described 
in the next section. 


Method of Solution 

We are dealing here essentially with a set of non-lin- 
ear coupled parabolic partial differential equations. A 

. . 7 

numerical scheme used previously by Sharma and Singnano 
for a similar system of equations has also been adopted 
in this particular case. We will describe below only sortie 
of the new features involved in the present problem. 

The flow outside the boundary layer is governed by a 
nonlinear ordinary differential equation, Eq. (32), and 

Q 

was solved by using fourth-order Runge-Kutta method. The 

initial condition, u g (x=0) , was determined by assuming one 

dimensional isentropic flow in the nozzle, and by using a 

value of Y„ which resulted in giving the pressure at the 
h 2 

leading edge of the specimen (p Q ) equal to that obtained 


by extrapolating the experimentally observed values. y 


H- 


obtained in this manner was always within l^f/o to a simple 

average of the -values corresponding to the chamber con- 

H 2 

ditions and to those commonly existing at the downstream 


end of the specimen. A second constant average value of 
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Y , which corresponded more closely to the temperature 
n 2 

variations in the domain of the specimen, was used in Eq. 

(32) , and was not changed from one test condition to an- 
other. These choices of Ytt led to the pressure profiles 

H 2 

in the domain of its specimen within 1 % of the more elab- 
orate two-dimensional, non-equilibrium calculations due 
9 — 

to Hamann. Thus, u g (x) was obtained along the specimen 

at stations separated by a uniform step size of Dx = 0.01. 

The other flow variables like n , p , T , etc. were then 

e e e 

obtained by using Eqs. (35) and (36) . 

The initial conditions for the governing boundary 
layer equations, Eqs. (20) to (23), were assumed to be given 
by a locally similar solution of these equations at a small 
enough value of s. A choice of = 0.01 was used and the 
results were found to be insensitive to any reasonable re- 
ductions in this value. The continuity equation, Eq. (20), 
having a first order derivative in r| was directly integrated 

g 

by using extended Simpson's rule. 


Note that the external flow field is determined at 


uniform intervals in x so that the corresponding values of 
s are no longer uniform, [see Eq. (15)]. Therefore, in 
order to obtain the difference-differential forms of Eqs. 

(20) to (23), we replace the s-derivat ives by the difference 

g 

expressions with the help of general Lagrange formula. 
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That is, 


au 

a s 


(k) 
(i/ j) 


U . . - U , . , .v 

(1, i) (1-1, i) 

s . - s . , 

1 l-l 


(42) 


for two-point approximation 


and 


/au 

\ds 


x (k) 

( i » j ) 


( Si - Si - 2 ) + (k) 

(s i- s i-2> 'V 8 !-! 1 (i ' j) 


(s i- s i-2> 


(s i- s i-i> (s i-r s i-2> 


U 


(i-l/ j) 


+ 




< s i- s i-2> <s i-l- s i-2> 


u 


( i-2 , j ) 


(43) 


for three-point approximation 


Two-point approximation, Eq. (42), will be used only once 
in moving the first step away from the locally similar solu- 
tion. The terms like -rr lc 

au \ 

c(|-U) + ("|^)('|^)}' before expressing them in the finite dif- 
ference approximations. All ri-derivat ives were replaced by 
central difference formulae. The various transport param- 
eters like C, Pr, Sc, etc. were determined from the results 
of the previous iteration. The linearization and iteration 

procedure was done in the same manner as described by Sharma 

. _ . . 7 

and Singnano. 


au\ 
a v 


are expanded into 


The step size in ri-direction was kept fixed at 0.05 
while that in x-direction varied from 0.05 for the first 
two stations to 0.10 till the end of the specimen. The 
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external flow conditions were found to be adequately satis- 
fied at T1 max = 7.00. The iteration cycle was stopped when 
a prescribed convergence of 1 in 10 4 was achieved for tem- 
perature values of successive iterations. 

Results and Discussion 

The detailed measurements in each test run of Refer- 
ence 2 consisted of the erosion rates, and the pressures 
at four locations along the specimen. In the present cal- 
culations, we obtain, in addition to the erosion rates and 
pressures, the detailed profiles of the x-component of ve- 
locity, temperature, and the species mass fractions as 
well as the rate of heat transfer at the specimen surface. 
Since the surface temperature of the specimen is not known 
a priori, a couple of trial runs are necessary to obtain 
a value of surface temperature which gives erosion rates 
close to the experimental points. 

The magnitudes of the various physical quantities in- 
volved in the computation are listed next. r fc = 0.3175 cm, 
a = 7°, r = 0.5144 cm, L = 3.95 cm, 7 = 1.385, W„ = 2.016 

^ Hn 

-1 -**1 A -1 
gm mole , W = 12.011 gm mole , and W = 16.043 gm mole . 

c CJrl . 

^ —2 _ i 

In order to transform the mass burning rate (gm cm sec ) 
into the measured erosion rate (cm sec -1 ) , we made use of 
the density of the AGCarb-101, that is, p c = 1.45 gm cm -3 , 
as determined at the Aerojet Nuclear Systems Company 10 . 

The results of the solution of Eqs. (20 to (23) for 
the stagnation conditions of test run No. W-110 of Reference 2 
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are plotted in Figs. 2 to 6 when = 2188°K. The agree- 
ment between the calculated and the measured pressure (see 
Fig. 2) is quite good, and the differences at various lo- 
cations ramain within a few percent. Note that an accurate 
determination of p is very important since the erosion rate 
is directly proportional to p [see Eq. (38)]. The general 
features of the velocity (Fig. 3) and the temperature (Fig. 4) 
profiles are similar to those commonly observed in other 
boundary layer flows. A plot of q w [ = (ST/By) w ] vs x in 
Fig. 5 shows that the heat is always being transferred from 
the gaseous phase to the specimen. Note that q^ is depen- 
dent on a choice for the initial value of x, and that the 
effect is insignificant towards the downstream portion of 
the specimen. A choice of it = 0.01 was found quite adequate 
for the present computations because it had avoided the sin- 
gularity at x = 0 as well as was reasonably small enough not to 
give results in appreciable error. The erosion rates were, 
however, insensitive to the value of x^ as long as we dealt 
with constant surface temperature conditions. Fig. 6 con- 
tains the comparison between the calculated erosion rates, 
and the experimental points. In view of some complexities 
involved in the experimental setup as well as the simpli- 
fications made in the present model, the agreement between 
the two results is considered quite satisfactory. It must 
be pointed out here that this comparison varies from one 
test run condition to another without any discernible 
pattern. Furthermore, since the experimental points for 



different test runs having similar stagnation conditions 
do not indicate similar functional dependences, there is 
a possibility of some errors in the measurements which 
have not been estimated. 

The amount of hydrogen gas consumed during the erosion 

process was found to be extremely small; Y„ being gener- 

H 2 / w 

ally of the order of 0.999. Thus, no detailed profiles for 
the species mass fractions were plotted. More importantly, 
the erosion rates can thus be calculated directly from Eq. 
(39) by setting Yjj w equal to unity, and by using the val- 
ues of p^ obtained from a solution of equations governing 
the external flow. Of course, the boundary layer equations 
must be solved if the heat transfer rates are to be deter- 
mined as well as if the surface temperature is to be deter- 
mined by solving the coupled heat transfer problem within 
the spec imen . 

The general features of the detailed profiles com- 
puted for other test conditions are very similar to those 
of the previous case except that q becomes negative for 
test conditions in which the (assumed) surface temperature 
is higher than the stagnation value. Some of the computed 
erosion rates for different test conditions are shown in 
Fig. 7. The agreement between the calculated curves and 
the experimental points is thus extremely good in tests 
like W-lll, W-112, and W-115, while there is an appreciable 
error for test conditions of W-142. Further improvements 
in the theoretical model are possible as well as some 



realistic error estimates of the measured quantities are 
necessary. 


Conclusions 

The agreement between the theoretical curves and the 
experimental points is quite reasonable, considering the 
fact that the surface boundary condition for the tempera- 
ture (or the energy) equation in the actual experiments is 
far more complicated than the one used in the present cal- 
culations. The general trend of the curves is also simi- 
lar to those of the experimental points. The surface re- 
action model, therefore, would appear adequate in the pre- 
diction of erosion rates for specimen surface temperatures 
less than 3000° K. Further improvements in the model can 
be achieved by including the dissociation and the turbu- 
lent flow effects as well as the coupled heat transfer 
problem within the specimen. In addition, for a precise 
agreement between theory and experiment, experiments with 
simplified energy boundary conditions, like constant sur- 
face temperature, must be performed. 

References 

'''Altseimer , J. H., Mader, G. F., and Stewart, J. J. , 
"Operating Characteristics and Requirements for the NERVA 
Flight Engine," AIAA Paper No. 70-676, June 1970, pp. 1-9. 

2 

Inouye, H. H., Claassen, L. B., and Shurley, L. A., 
"Final Test Report: AGCarb-101 Corrosion/Erosion Evaluation 

Program," 7750:R0834, April 1970, Aerojet Nuclear Systems 
Company, Sacramento, Calif. 



Reaction of Graphite 


^Clarke, J. T, and Fox, B. R., 

Filaments with Hydrogen above 2000°K, "The Journal of 
Chem. Phys., Vol. 46, No. 3, Feb. 1967, pp. 827-836. 

^Lees, L. , "Laminar Heat Transfer over Blunt-Nosed 
Bodies at Hypersonic Flight Speeds," Jet Propulsion, Vol. 

26, No. 4, April 1956, pp. 259-269. 

^McBride, B. J., Heimel, S., Ehlers, J. G., and Gorden, S. 
"Thermodynamic Properties to 6000° K for 210 Substances In- 
volving the First 18 Elements," NASA SP-3001, March 1963. 

6penner, S. S., "Chemistry Problems in Jet Propulsion," 
MacMillan, New York, 1957, pp. 245-250. 

^Sharma, 0. P. and Sirignano, W. A., "On the Ignition 
of a Pre-Mixed Fuel by a Hot Projectile," Comb. Sc. and 
Tech., Vol. 1, No. 6, May 1970, pp. 481-494. 

Q 

Abramowitz , . M. and Stegun, I. A., ed. , "Handbook of 
Mathematical Functions with Formulas, Graphs, and Mathe- 
matical Tables," NBS AMS 55, March 1965, p. 896, p. 886, 
pp. 882-884. 

9 Hamann, R. J., private communication, Guggenheim 
Laboratories, Princeton University, Princeton, New Jersey. 

l°David, H. 0., Kotfila, R. J. , and Schleicher, W. J., 

"The Fabrication and Properties of the Fibrous Reinforced 
Graphite Composite, AGCarb-101, " RN-S-0549, March 1970, 
Aerojet Nuclear Systems Company, Sacramento, California. 



19 


Nomenc lature 

a known constant [see Eq. (39)] (atm) 
specific heat at constant pressure of species i 
(cal gm K 


(PU / P e U e ) 

2 -1 

binary diffusion coefficient (cm sec ) 

a known constant [see Eq. (39)] (cal mole 1 ) 

T r 

specific enthalpy of species i [J c ^ dT + h^ 

ip 1? > 

r 


] (cal gm 1 ) 


thermal conductivity (cal cm ^sec ^°K ^) 
dissociation constant [see Eq. (39)] (atm) 

_2 

known constant [see Eq. (38) ] (moles of carbon cm sec 

(atm) -3/ " 2 ) 


length of the graphite specimen (cm) 
pressure (atm) 

Prandtl number [= |j (S c .Y.)/k] 

i 

heat flux [= k (?>T/dy ) ] (cal cm 2 sec ^) 
perpendicular distance of the specimen surface from 
the axis of symmetry (cm) 
radius of the nozzle at the leading edge of the 
graphite specimen (cm) 
radius of the nozzle throat (cm) 
universal gas constant (cal mole ^ °K ^ ) 

modified coordinate in x-direction [see Eq. (15)] 

2 -2 
(gm sec ) 

Schmidt number [= (p./ pD ) ] 

temperature (°K) 
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u = x-component of velocity (cm sec”^) 

U = non-dimensional velocity [see Eq. (17)] 

v = y-component of velocity (cm sec 

V = modified transverse velocity [see Eq. (19)] 

W i = molecular weight of species i (gm mole -1 ) 

x = distance along the specimen surface (cm) 

x = non-dimensional distance [=(x/L)] 

y = distance perpendicular to the specimen surface (cm) 

= mass fraction of species i 

a = half angle of the divergent part of the nozzle [= 7°] 

p = (2s/u e ) (du e /ds) 

Y = ratio of the specific heats at constant pressure 

and at constant volume 
T = (2s/T g ) (d T 0 /ds) 

T| = transverse coordinate [see Eq. (16)] 

9 = dimensionless temperature [see Eq. (18)] 

^ = coefficient of viscosity (gm cm -1 sec -1 ) 

_3 

p = density (gm cm ) 

• -2 -1 
uo ^ = reaction rate of species i (gm cm sec ) 

Subscripts 

c = carbon 

e = external flow (y=«) 

i = species (H 2 and CH 4 ) or a mesh point in s-direction 

j = a mesh point in ri-direction 

= leading edge of the specimen (x = 0) 


o 



r 


reference state 


s = specimen 

w = specimen surface (y = 0) 

Superscripts 

= non-dimensional or average 
k = iteration number 
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